#include  "solver.h"
#include  "transport.h"
#include  "burgers.h"
#include  "euler.h"

/*--------------------------------------------------------*/

int main(int argc, char** argv)
{
  if(argc != 6)
  {
    std::cerr << "usage: " << argv[0] << " problem nnodes degre dt finaltime\n";
    exit(1);
  }

  ProblemInterface* problem = NULL;
  std::string problemname(argv[1]);
  if(problemname == "Transport")
  {
    problem = new Transport;
  }
  else if(problemname == "Burgers")
  {
    problem = new Burgers;
  }
  else if(problemname == "Euler")
  {
    problem = new Euler;
  }
  else
  {
    std::cerr << "*** ERROR in " << argv[0] << " : unknown problem \"" << problemname << "\"\n";
    exit(1);
  }
  int nnodes = atoi(argv[2]);
  int degre = atoi(argv[3]);
  double dt = atof(argv[4]);
  double finaltime = atof(argv[5]);
    
  Mesh mesh(nnodes);
  int ncells = mesh.ncells();

  int ncomp = problem->ncomp();
  SystemVector u(ncomp, degre), uold(ncomp, degre), f(ncomp, degre);

  u.reSize(ncells);
  uold.reSize(ncells);
  f.reSize(ncells);

  Solver S(degre);

  std::string filename("solution");

  int niter = (int) (finaltime/dt);
  double time = 0.0;

  S.initialCondition(u, problem, mesh);

  // std::cerr << "u " << u << "\n";

  S.writeSolution(filename, time, u, problem, mesh);

  for(int iter = 0; iter < niter; iter++)
  {
    uold = u;
    f.zero();
    S.computeResidual(f, uold, problem, mesh);

    for(int icomp = 0; icomp < u.ncomp(); icomp++)
    {
      for(int icell = 0; icell < ncells; icell++)
      {
        double x1, x2;
        mesh.nodesOfCell(x1, x2, icell);
        double dx = x2-x1;
        for(int i = 0; i < degre+1; i++)
        {
          u(icomp) (icell, i) -= ( dt/dx )*f(icomp) (icell, i);
        }
      }
    }
    time += dt;
    S.writeSolution(filename, time, u, problem, mesh);
  }
  S.writeSolution(filename, u, problem, mesh);

  delete problem;
  return 0;
}

